###############################################################################################################
#The compilation and creation of the data for "How Parliamentarism Developed in Western Europe" (ParlDev data)#
###############################################################################################################


setwd("C:/Users/simda81/OneDrive - Linköpings universitet/Parliamentary Development/Simon avhandling/BJPS/Replication files and data")


#library(ggplot2)
library(plyr)
library(tidyverse)
library(viridis)
library(hrbrthemes)
library(boot)
library(car)
library(bsts)
library(haven)
library(here)

####List of variables produced

### mean04=parliamentarism, main estimation
### mean03=parliamentarism, decreased uncertainty paramenter
### mean05=parliamentarism, increased uncertainty parameter
### meanconst = estimation adjusted for major constitutional reforms or regime changes
### meanhosstart = estimation assuming all countries were personal monarchies at the start,
# i.e., that only the head of state is expected to terminate the government
### meannoaftt=estimation that disregards failed termination attempts (no other actor fails to terminate)
### meannoother=estimation that disregards "other" events, i.e., those in which neither parliament
# nor the head of state is involved.

#####Individual country graphs

###Austria

austrias<-read.csv2(file="Austriashort.csv", fileEncoding = "ISO-8859-1")

austrianew <- data.frame(year=austrias$year[1]:austrias$year[nrow(austrias)])

austria <- full_join(austrias,austrianew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

austria <- austria %>%  add_row(year = full_seq(c(1861, 1864), 1), .before = 1) 

austria$event <- with(austria, ifelse(is.na(event), 0, event))

austria$eventnoaftt<-with(austria, ifelse(year == 1879 & government == "Taaffe" | year == 1915, NA, noncevent))

austria$noother<-with(austria, ifelse(noncevent == 1 & pevent == 0, NA, pevent))

###SdPrior = 0.3

austriass03<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.3, 
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_complete03<- bsts(austria$noncevent,
                             state.specification=austriass03,
                             niter=10000,
                             family="logit",
                             seed=123456,
                             save.full.state="TRUE",
)

modelaustria_complete03$timestamp.info$timestamps <- austria$year 
modelaustria_complete03$timestamp.info$regular.timestamps <- austria$year 


###SdPrior = 0.5

austriass05<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.5, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_complete05<- bsts(austria$noncevent,
                              state.specification=austriass05,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelaustria_complete05$timestamp.info$timestamps <- austria$year 
modelaustria_complete05$timestamp.info$regular.timestamps <- austria$year 

###SdPrior = 0.4


austriass04<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.4, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_complete04<- bsts(austria$noncevent,
                              state.specification=austriass04,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelaustria_complete04$timestamp.info$timestamps <- austria$year 
modelaustria_complete04$timestamp.info$regular.timestamps <- austria$year 

### With SdPrior=0.6 (for constitutional changes)

austriassconst<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.6, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_completeconst<- bsts(austria$noncevent,
                              state.specification=austriassconst,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelaustria_completeconst$timestamp.info$timestamps <- austria$year 
modelaustria_completeconst$timestamp.info$regular.timestamps <- austria$year 

#Hosstart
austriassJ<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.4, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_completeJ<- bsts(austria$noncevent,
                              state.specification=austriassJ,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelaustria_completeJ$timestamp.info$timestamps <- austria$year 
modelaustria_completeJ$timestamp.info$regular.timestamps <- austria$year 

###noaftt

austriassnoaftt<-AddLocalLevel(list(),  
                               sigma.prior=SdPrior(sigma.guess=.4, 
                                                   sample.size=10,
                                                   upper.limit=Inf),
                               initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_completenoaftt<- bsts(austria$eventnoaftt,
                                   state.specification=austriassnoaftt,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
                                   save.full.state="TRUE",
)

modelaustria_completenoaftt$timestamp.info$timestamps <- austria$year 
modelaustria_completenoaftt$timestamp.info$regular.timestamps <- austria$year 

###noother

austriassnoother<-AddLocalLevel(list(),  
                               sigma.prior=SdPrior(sigma.guess=.4, 
                                                   sample.size=10,
                                                   upper.limit=Inf),
                               initial.state.prior=NormalPrior(-7, 1)) 


modelaustria_completenoother<- bsts(austria$noother,
                                   state.specification=austriassnoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
                                   save.full.state="TRUE",
)

modelaustria_completenoother$timestamp.info$timestamps <- austria$year 
modelaustria_completenoother$timestamp.info$regular.timestamps <- austria$year 


modelaustria_complete04$state.contributions[1:10000,1,which(austria$year==1933):which(austria$year==1944)] <- NA

pdf("austria.pdf")
plot(modelaustria_complete04,  ylim=c(0,1), scale="mean",
     #main="Austria",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

png("austria.png", width=17,height=17,units="cm", res=600)
plot(modelaustria_complete04,  ylim=c(0,1), scale="mean",
     #main="Austria",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()




###Belgium



belgiums<-read.csv2(file="Belgiumshort.csv", fileEncoding = "ISO-8859-1")


belgiumnew <- data.frame(year=belgiums$year[1]:belgiums$year[nrow(belgiums)])

belgium <- full_join(belgiums,belgiumnew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))


belgium <- add_row(belgium, country = "Belgium", year = 2021, event = 0) #Add years for Belgium 

belgium$event <- with(belgium, ifelse(is.na(event), 0, event)) #To get the events right

belgium$eventnoaftt<-with(belgium, ifelse(year == 1940 | year == 1960, NA, noncevent))


belgium$noother<-with(belgium, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


#With SdPrior = 0.3

belgiumss03<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.3, 
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-2, 1)) 


modelbelgium_complete03<- bsts(belgium$noncevent,
                             state.specification=belgiumss03,
                             niter=10000,
                             family="logit",
                             seed=123456,
                             save.full.state="TRUE",
)

modelbelgium_complete03$timestamp.info$timestamps <- belgium$year 
modelbelgium_complete03$timestamp.info$regular.timestamps <- belgium$year 

### With SdPrior=0.4

belgium04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=0.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-2, 1))  


modelbelgium_complete04<- bsts(belgium$noncevent,
                              state.specification=belgium04,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelbelgium_complete04$timestamp.info$timestamps <- belgium$year 
modelbelgium_complete04$timestamp.info$regular.timestamps <- belgium$year 

### With SdPrior=0.5

belgium05<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=0.5, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-2, 1))  


modelbelgium_complete05<- bsts(belgium$noncevent,
                              state.specification=belgium05,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelbelgium_complete05$timestamp.info$timestamps <- belgium$year 
modelbelgium_complete05$timestamp.info$regular.timestamps <- belgium$year 

#Hosstart

belgiumJ<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=0.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-7, 1)) 


modelbelgium_completeJ<- bsts(belgium$noncevent,
                              state.specification=belgiumJ,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              save.full.state="TRUE",
)

modelbelgium_completeJ$timestamp.info$timestamps <- belgium$year 
modelbelgium_completeJ$timestamp.info$regular.timestamps <- belgium$year 



### Noaftt

belgiumnoaftt<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=0.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(-2, 1)) 


modelbelgium_completenoaftt<- bsts(belgium$eventnoaftt,
                                   state.specification=belgiumnoaftt,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
                                   save.full.state="TRUE",
)

modelbelgium_completenoaftt$timestamp.info$timestamps <- belgium$year 
modelbelgium_completenoaftt$timestamp.info$regular.timestamps <- belgium$year 

#noother

belgiumnoother<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=0.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(-2, 1)) 


modelbelgium_completenoother<- bsts(belgium$noother,
                                   state.specification=belgiumnoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
                                   save.full.state="TRUE",
)

modelbelgium_completenoother$timestamp.info$timestamps <- belgium$year 
modelbelgium_completenoother$timestamp.info$regular.timestamps <- belgium$year 


###Graphs

png("belgium.png", width=17,height=17,units="cm", res=600)
plot(modelbelgium_complete04,  ylim=c(0,1), scale="mean",
     #     main="Belgium",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()


###Britain

britains<-read.csv2(file="Britainshort.csv", fileEncoding = "ISO-8859-1")

britainnew <- data.frame(year=britains$year[1]:britains$year[nrow(britains)])

britain <- full_join(britains,britainnew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))


britain <- add_row(britain, country = "Britain", year = 2020, event = 0) #Add years for Britain
britain <- add_row(britain, country = "Britain", year = 2021, event = 0) #Add years for Britain


britain$event <- with(britain, ifelse(is.na(event), 0, event))

britain$eventnoaftt<-with(britain, ifelse(year == 1784, NA, noncevent))

britain$noother<-with(britain, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###With SdPrior = 0.3

britainss03<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.3, 
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-1, 1)) 

modelbritain_complete03<- bsts(britain$noncevent,
                             state.specification=britainss03,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelbritain_complete03$timestamp.info$timestamps <- britain$year 
modelbritain_complete03$timestamp.info$regular.timestamps <- britain$year 

###With SdPrior=0.4


britain04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-1, 1)) 

modelbritain_complete04<- bsts(britain$noncevent,
                              state.specification=britain04,
                              niter=10000,
                              family="logit",
                              seed=123456,
)


modelbritain_complete04$timestamp.info$timestamps <- britain$year 
modelbritain_complete04$timestamp.info$regular.timestamps <- britain$year 


###With SdPrior=0.5


britain05<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.5, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-1, 1)) 

modelbritain_complete05<- bsts(britain$noncevent,
                              state.specification=britain05,
                              niter=10000,
                              family="logit",
                              seed=123456,
)


modelbritain_complete05$timestamp.info$timestamps <- britain$year 
modelbritain_complete05$timestamp.info$regular.timestamps <- britain$year 



#Hosstart

britainJ<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modelbritain_completeJ<- bsts(britain$noncevent,
                              state.specification=britainJ,
                              niter=10000,
                              family="logit",
                              seed=123456,
)


modelbritain_completeJ$timestamp.info$timestamps <- britain$year 
modelbritain_completeJ$timestamp.info$regular.timestamps <- britain$year 


###Noaftt


britainnoaftt<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(-1, 1)) 

modelbritain_completenoaftt<- bsts(britain$eventnoaftt,
                                   state.specification=britainnoaftt,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)


modelbritain_completenoaftt$timestamp.info$timestamps <- britain$year 
modelbritain_completenoaftt$timestamp.info$regular.timestamps <- britain$year 



###Noother


britainnoother<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(-1, 1)) 

modelbritain_completenoother<- bsts(britain$noother,
                                   state.specification=britainnoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)


modelbritain_completenoother$timestamp.info$timestamps <- britain$year 
modelbritain_completenoother$timestamp.info$regular.timestamps <- britain$year 


###Country graphs


png("britain.png", width=17,height=17,units="cm", res=600)
plot(modelbritain_complete04,  ylim=c(0,1), scale="mean",
     #main="Britain",# sigma.guess=0.4, -2500, 50", 
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()


###Denmark


denmarks<-read.csv2(file="Denmarkshort.csv", fileEncoding = "ISO-8859-1")


denmarknew <- data.frame(year=denmarks$year[1]:denmarks$year[nrow(denmarks)])

denmark <- full_join(denmarks,denmarknew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

denmark <- denmark %>%  add_row(year = full_seq(c(1849, 1851), 1), .before = 1) #DATA SENSITIVE


denmark <- add_row(denmark, country = "Denmark", year = 2020, event = 0) #Add years for Denmark
denmark <- add_row(denmark, country = "Denmark", year = 2021, event = 0) #Add years for Denmark


denmark$event <- with(denmark, ifelse(is.na(event), 0, event))

denmark$eventnoaftt <- with(denmark, ifelse(year == 1854 & what == "An address of no confidence from both houses against the government to the king is unheeded" |
                                              year == 1873 | year == 1883 | year == 1885, NA, noncevent))

denmark$noother<-with(denmark, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###With SdPrior = 0.3

denmarkss03<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.3,
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(0, 1)) 

modeldenmark_complete03<- bsts(denmark$noncevent,
                             state.specification=denmarkss03,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modeldenmark_complete03$timestamp.info$timestamps <- denmark$year 
modeldenmark_complete03$timestamp.info$regular.timestamps <- denmark$year 

#With SdPrior=0.4

denmark04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4,
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(0, 1)) 

modeldenmark_complete04<- bsts(denmark$noncevent,
                              state.specification=denmark04,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modeldenmark_complete04$timestamp.info$timestamps <- denmark$year 
modeldenmark_complete04$timestamp.info$regular.timestamps <- denmark$year 

#With SdPrior=0.5

denmark05<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.5,
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(0, 1)) 

modeldenmark_complete05<- bsts(denmark$noncevent,
                              state.specification=denmark05,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modeldenmark_complete05$timestamp.info$timestamps <- denmark$year 
modeldenmark_complete05$timestamp.info$regular.timestamps <- denmark$year 

#Hosstart
denmarkJ<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4,
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modeldenmark_completeJ<- bsts(denmark$noncevent,
                              state.specification=denmarkJ,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modeldenmark_completeJ$timestamp.info$timestamps <- denmark$year 
modeldenmark_completeJ$timestamp.info$regular.timestamps <- denmark$year 



#Noaftt

denmarknoaftt<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4,
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(0, 1)) 

modeldenmark_completenoaftt<- bsts(denmark$eventnoaftt,
                                   state.specification=denmarknoaftt,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)

modeldenmark_completenoaftt$timestamp.info$timestamps <- denmark$year 
modeldenmark_completenoaftt$timestamp.info$regular.timestamps <- denmark$year 

#Noother

denmarknoother<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4,
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(0, 1)) 

modeldenmark_completenoother<- bsts(denmark$noother,
                                   state.specification=denmarknoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)

modeldenmark_completenoother$timestamp.info$timestamps <- denmark$year 
modeldenmark_completenoother$timestamp.info$regular.timestamps <- denmark$year 


#Country graphs

png("Denmark.png", width=17,height=17,units="cm", res=600)
plot(modeldenmark_complete03,  ylim=c(0,1), scale="mean",
     #main="Denmark", #, sigma.guess=0.4, -100, 50", 
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

###Finland



finlands<-read.csv2(file="Finlandshort.csv", fileEncoding = "ISO-8859-1")

finlandnew <- data.frame(year=finlands$year[1]:finlands$year[nrow(finlands)])

finland <- full_join(finlands,finlandnew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

finland <- finland %>%  add_row(year = full_seq(c(1863, 1881), 1), .before = 1) #DATA SENSITIVE

finland <- add_row(finland, country = "Finland", year = 2020, event = 0) #Add years for Finland
finland <- add_row(finland, country = "Finland", year = 2021, event = 0) #Add years for Finland


finland$event <- with(finland, ifelse(is.na(event), 0, event))

finland$eventnoaftt<-finland$noncevent

finland$noother<-with(finland, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###SdPrior = 0.3

finlandss03<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=0.3, 
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-5, 1)) 


modelfinland_complete03<- bsts(finland$noncevent,
                             state.specification=finlandss03,
                             niter=10000,
                             family="logit",
                             seed=123456,
                             #save.full.state="TRUE",
)

modelfinland_complete03$timestamp.info$timestamps <- finland$year 
modelfinland_complete03$timestamp.info$regular.timestamps <- finland$year 

###Finland 0.4

finlandss04<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=0.4, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-5, 1)) 


modelfinland_complete04<- bsts(finland$noncevent,
                              state.specification=finlandss04,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              #save.full.state="TRUE",
)

modelfinland_complete04$timestamp.info$timestamps <- finland$year 
modelfinland_complete04$timestamp.info$regular.timestamps <- finland$year 

###sdPrior = 0.5

finland05<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=0.5, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 


modelfinland_complete05<- bsts(finland$noncevent,
                              state.specification=finland05,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              #save.full.state="TRUE",
)

modelfinland_complete05$timestamp.info$timestamps <- finland$year 
modelfinland_complete05$timestamp.info$regular.timestamps <- finland$year 

#Hosstart
finlandssJ<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=0.4, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-5, 1)) 


modelfinland_completeJ<- bsts(finland$noncevent,
                              state.specification=finlandssJ,
                              niter=10000,
                              family="logit",
                              seed=123456,
                              #save.full.state="TRUE",
)

modelfinland_completeJ$timestamp.info$timestamps <- finland$year 
modelfinland_completeJ$timestamp.info$regular.timestamps <- finland$year 

###Finland noother

finlandssnoother<-AddLocalLevel(list(),  
                           sigma.prior=SdPrior(sigma.guess=0.4, 
                                               sample.size=10,
                                               upper.limit=Inf),
                           initial.state.prior=NormalPrior(-5, 1)) 


modelfinland_completenoother<- bsts(finland$noother,
                               state.specification=finlandssnoother,
                               niter=10000,
                               family="logit",
                               seed=123456,
                               #save.full.state="TRUE",
)

modelfinland_completenoother$timestamp.info$timestamps <- finland$year 
modelfinland_completenoother$timestamp.info$regular.timestamps <- finland$year 

#Country graphs

png("finland.png", width=17,height=17,units="cm", res=600)
plot(modelfinland_complete04,  ylim=c(0,1), scale="mean",
     #main="Finland",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

###France

frances<-read.csv2(file="Franceshort.csv", fileEncoding = "ISO-8859-1")


francenew <- data.frame(year=frances$year[1]:frances$year[nrow(frances)])

france <- full_join(frances,francenew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))


france <- add_row(france, country = "France", year = 2021, event = 0) #Add years for France

france$event <- with(france, ifelse(is.na(event), 0, event))

france$eventnoaftt<-with(france, ifelse(government == "Rochebouet" & nonpevent == 1, NA, noncevent))

france$noother<-with(france, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###SdPrior = 0.3

francess03<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.3, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete03<- bsts(france$noncevent,
                            state.specification=francess03,
                            niter=10000,
                            family="logit",
                            seed=123456,
)

modelfrance_complete03$timestamp.info$timestamps <- france$year 
modelfrance_complete03$timestamp.info$regular.timestamps <- france$year 

#With SdPrior=0.4

france04<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete04<- bsts(france$noncevent,
                             state.specification=france04,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelfrance_complete04$timestamp.info$timestamps <- france$year 
modelfrance_complete04$timestamp.info$regular.timestamps <- france$year 

#With SdPrior=0.5

francess05<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.5, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete05<- bsts(france$noncevent,
                             state.specification=francess05,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelfrance_complete05$timestamp.info$timestamps <- france$year 
modelfrance_complete05$timestamp.info$regular.timestamps <- france$year 

#France with SdPrior = 0.9

franceconst<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.9, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completeconst<- bsts(france$noncevent,
                             state.specification=franceconst,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelfrance_completeconst$timestamp.info$timestamps <- france$year 
modelfrance_completeconst$timestamp.info$regular.timestamps <- france$year 

#Hosstart
franceJ<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-7, 1)) 

modelfrance_completeJ<- bsts(france$noncevent,
                             state.specification=franceJ,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelfrance_completeJ$timestamp.info$timestamps <- france$year 
modelfrance_completeJ$timestamp.info$regular.timestamps <- france$year 


#noaftt

francenoaftt<-AddLocalLevel(list(),  
                            sigma.prior=SdPrior(sigma.guess=.4, 
                                                sample.size=10,
                                                upper.limit=Inf),
                            initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completenoaftt<- bsts(france$eventnoaftt,
                                  state.specification=francenoaftt,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)

modelfrance_completenoaftt$timestamp.info$timestamps <- france$year 
modelfrance_completenoaftt$timestamp.info$regular.timestamps <- france$year 


#noother

francenoother<-AddLocalLevel(list(),  
                            sigma.prior=SdPrior(sigma.guess=.4, 
                                                sample.size=10,
                                                upper.limit=Inf),
                            initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completenoother<- bsts(france$noother,
                                  state.specification=francenoother,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)

modelfrance_completenoother$timestamp.info$timestamps <- france$year 
modelfrance_completenoother$timestamp.info$regular.timestamps <- france$year 


### Country graphs



modelfrance_complete04$state.contributions[1:10000,1,which(france$year==1940):which(france$year==1944)] <- NA

png("france.png", width=17,height=17,units="cm", res=600)
plot(modelfrance_complete04,  ylim=c(0,1), scale="mean",
     #main="France",#, sigma.guess=0.4, -3500, 50", 
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

###Germany

germanys<-read.csv2(file="Germanyshort.csv", fileEncoding = "ISO-8859-1")


germanynew <- data.frame(year=germanys$year[1]:germanys$year[nrow(germanys)]) 

germany <- full_join(germanys,germanynew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

germany <- germany %>%  add_row(year = full_seq(c(1871, 1889), 1), .before = 1) #DATA SENSITIVE


germany$event <- with(germany, ifelse(is.na(event), 0, event))

germany$eventnoaftt <- with(germany, ifelse(year == 1913, NA, noncevent))

germany$noother<-with(germany, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###SdPrior = 0.5

germanyss05<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.5,
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_complete05<- bsts(germany$noncevent,
                             state.specification=germanyss05,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelgermany_complete05$timestamp.info$timestamps <- germany$year 
modelgermany_complete05$timestamp.info$regular.timestamps <- germany$year 

#With SdPrior = 0.3

germany03<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.3, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_complete03<- bsts(germany$noncevent,
                              state.specification=germany03,
                              niter=10000,
                              family="logit",
                              seed=123456
)


modelgermany_complete03$timestamp.info$timestamps <- germany$year 
modelgermany_complete03$timestamp.info$regular.timestamps <- germany$year 

#With SdPrior = 0.4

germany04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_complete04<- bsts(germany$noncevent,
                              state.specification=germany04,
                              niter=10000,
                              family="logit",
                              seed=123456
)


modelgermany_complete04$timestamp.info$timestamps <- germany$year 
modelgermany_complete04$timestamp.info$regular.timestamps <- germany$year 

# with SdPrior = 0.6 (for constitutional changes)


germanyconst<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.6, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_completeconst<- bsts(germany$noncevent,
                              state.specification=germanyconst,
                              niter=10000,
                              family="logit",
                              seed=123456
)


modelgermany_completeconst$timestamp.info$timestamps <- germany$year 
modelgermany_completeconst$timestamp.info$regular.timestamps <- germany$year 

#Hosstart (same as baseline)

germanyJ<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_completeJ<- bsts(germany$noncevent,
                              state.specification=germanyJ,
                              niter=10000,
                              family="logit",
                              seed=123456
)


modelgermany_completeJ$timestamp.info$timestamps <- germany$year 
modelgermany_completeJ$timestamp.info$regular.timestamps <- germany$year 

#Noaftt

germanyssnoaftt<-AddLocalLevel(list(),  
                               sigma.prior=SdPrior(sigma.guess=.4, 
                                                   sample.size=10,
                                                   upper.limit=Inf),
                               initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_completenoaftt<- bsts(germany$eventnoaftt,
                                   state.specification=germanyssnoaftt,
                                   niter=10000,
                                   family="logit",
                                   seed=123456
)


modelgermany_completenoaftt$timestamp.info$timestamps <- germany$year 
modelgermany_completenoaftt$timestamp.info$regular.timestamps <- germany$year 

#Noother

germanyssnoother<-AddLocalLevel(list(),  
                               sigma.prior=SdPrior(sigma.guess=.4, 
                                                   sample.size=10,
                                                   upper.limit=Inf),
                               initial.state.prior=NormalPrior(-5, 1)) 

modelgermany_completenoother<- bsts(germany$noother,
                                   state.specification=germanyssnoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456
)


modelgermany_completenoother$timestamp.info$timestamps <- germany$year 
modelgermany_completenoother$timestamp.info$regular.timestamps <- germany$year 

### Country graphs

modelgermany_complete04$state.contributions[1:10000,1,which(germany$year==1934):which(germany$year==1948)] <- NA

png("germany.png", width=17,height=17,units="cm", res=600)
plot(modelgermany_complete04,  ylim=c(0,1), scale="mean",
     #main="Germany",# sigma.guess=0.4", 
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()



###Italy



italys<-read.csv2(file="Italyshort.csv", fileEncoding = "ISO-8859-1")

italynew <- data.frame(year=italys$year[1]:italys$year[nrow(italys)]) 

italy <- full_join(italys,italynew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))


italy$event <- with(italy, ifelse(is.na(event), 0, event))

italy$eventnoaftt <- with(italy, ifelse(year == 1867 & government == "Menabrea", NA, noncevent))

italy$noother<-with(italy, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###With SdPrior = 0.4

italyss04<-AddLocalLevel(list(), 
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(1, 1)) 


modelitaly_complete04<- bsts(italy$noncevent,
                           state.specification=italyss04,
                           niter=10000,
                           family="logit",
                           seed=123456,
)


modelitaly_complete04$timestamp.info$timestamps <- italy$year 
modelitaly_complete04$timestamp.info$regular.timestamps <-italy$year 

# With SdPrior = 0.3

italy03<-AddLocalLevel(list(),  
                      sigma.prior=SdPrior(sigma.guess=.3, 
                                          sample.size=10,
                                          upper.limit=Inf),
                      initial.state.prior=NormalPrior(1, 1)) 


modelitaly_complete03<- bsts(italy$noncevent,
                            state.specification=italy03,
                            niter=10000,
                            family="logit",
                            seed=123456,
)


modelitaly_complete03$timestamp.info$timestamps <- italy$year 
modelitaly_complete03$timestamp.info$regular.timestamps <-italy$year 

#With SdPrior = 0.5

italy05<-AddLocalLevel(list(),  
                      sigma.prior=SdPrior(sigma.guess=.5, 
                                          sample.size=10,
                                          upper.limit=Inf),
                      initial.state.prior=NormalPrior(1, 1)) 


modelitaly_complete05<- bsts(italy$noncevent,
                            state.specification=italy05,
                            niter=10000,
                            family="logit",
                            seed=123456,
)


modelitaly_complete05$timestamp.info$timestamps <- italy$year 
modelitaly_complete05$timestamp.info$regular.timestamps <-italy$year 

#Hosstart

italyssJ<-AddLocalLevel(list(), 
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-6, 1)) 


modelitaly_completeJ<- bsts(italy$noncevent,
                            state.specification=italyssJ,
                            niter=10000,
                            family="logit",
                            seed=123456,
)


modelitaly_completeJ$timestamp.info$timestamps <- italy$year 
modelitaly_completeJ$timestamp.info$regular.timestamps <-italy$year 

#Noaftt


italyssnoaftt<-AddLocalLevel(list(), 
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(1, 1)) 


modelitaly_completenoaftt<- bsts(italy$eventnoaftt,
                                 state.specification=italyssnoaftt,
                                 niter=10000,
                                 family="logit",
                                 seed=123456,
)


modelitaly_completenoaftt$timestamp.info$timestamps <- italy$year 
modelitaly_completenoaftt$timestamp.info$regular.timestamps <-italy$year 

#Noother


italyssnoother<-AddLocalLevel(list(), 
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(1, 1)) 


modelitaly_completenoother<- bsts(italy$noother,
                                 state.specification=italyssnoother,
                                 niter=10000,
                                 family="logit",
                                 seed=123456,
)

modelitaly_completenoother$timestamp.info$timestamps <- italy$year 
modelitaly_completenoother$timestamp.info$regular.timestamps <-italy$year 

###Country graphs

modelitaly_complete03$state.contributions[1:10000,1,which(italy$year==1923):which(italy$year==1944)] <- NA

png("italy.png", width=17,height=17,units="cm", res=600)
plot(modelitaly_complete03,  ylim=c(0,1), scale="mean",
     #main="Italy",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

###Netherlands


netherlandss<-read.csv2(file="Netherlandsshort.csv", fileEncoding = "ISO-8859-1")

netherlandsnew <- data.frame(year=netherlandss$year[1]:netherlandss$year[nrow(netherlandss)]) 

netherlandsnew_join <- full_join(netherlandss,netherlandsnew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))


#The Netherlands before 1814 do not seem to meet the scope conditions, particularly not for any longer, continuous periods.
#The government seems to be a rotating parliamentary committee with a nominal head of government and state in one.
#The early nineteenth century, before French annexation are not meaningful to study from the general perspective applied.
netherlands<-subset.data.frame(netherlandsnew_join, netherlandsnew_join$year>1813) 

netherlands <- netherlands %>%  add_row(year = full_seq(c(2013, 2014), 1), .after = 204) #DATA SENSITIVE

netherlands <- netherlands %>%  add_row(year = 2016, .after = 207)

netherlands <- netherlands %>%  add_row(year = 2018, .after = 209)

netherlands <- netherlands %>%  add_row(year = full_seq(c(2020, 2021), 1), .after = 211)

netherlands$event <- with(netherlands, ifelse(is.na(event), 0, event))

netherlands$eventnoaftt <-netherlands$noncevent

netherlands$noother<-with(netherlands, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


###SdPrior = 0.4

netherlandsss04<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf), 
                             initial.state.prior=NormalPrior(-3, 1))


modelnetherlands_complete04<- bsts(netherlands$noncevent,
                                 state.specification=netherlandsss04,
                                 niter=10000,
                                 family="logit",
                                 seed=123456,
)


modelnetherlands_complete04$timestamp.info$timestamps <- netherlands$year 
modelnetherlands_complete04$timestamp.info$regular.timestamps <-netherlands$year 

#SdPrior = 0.3

netherlands03<-AddLocalLevel(list(),  
                            sigma.prior=SdPrior(sigma.guess=.3, 
                                                sample.size=10,
                                                upper.limit=Inf), 
                            initial.state.prior=NormalPrior(-3, 1))


modelnetherlands_complete03<- bsts(netherlands$noncevent,
                                  state.specification=netherlands03,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelnetherlands_complete03$timestamp.info$timestamps <- netherlands$year 
modelnetherlands_complete03$timestamp.info$regular.timestamps <-netherlands$year 

#SdPrior = 0.5

netherlands05<-AddLocalLevel(list(),  
                            sigma.prior=SdPrior(sigma.guess=.5, 
                                                sample.size=10,
                                                upper.limit=Inf), 
                            initial.state.prior=NormalPrior(-3, 1))


modelnetherlands_complete05<- bsts(netherlands$noncevent,
                                  state.specification=netherlands05,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelnetherlands_complete05$timestamp.info$timestamps <- netherlands$year 
modelnetherlands_complete05$timestamp.info$regular.timestamps <-netherlands$year 


#Hosstart
netherlandsssJ<-AddLocalLevel(list(),  
                              sigma.prior=SdPrior(sigma.guess=.4, 
                                                  sample.size=10,
                                                  upper.limit=Inf), 
                              initial.state.prior=NormalPrior(-3, 1))


modelnetherlands_completeJ<- bsts(netherlands$noncevent,
                                  state.specification=netherlandsssJ,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelnetherlands_completeJ$timestamp.info$timestamps <- netherlands$year 
modelnetherlands_completeJ$timestamp.info$regular.timestamps <-netherlands$year 

###noother

netherlandsssnoother<-AddLocalLevel(list(),  
                               sigma.prior=SdPrior(sigma.guess=.4, 
                                                   sample.size=10,
                                                   upper.limit=Inf), 
                               initial.state.prior=NormalPrior(-3, 1))


modelnetherlands_completenoother<- bsts(netherlands$noother,
                                   state.specification=netherlandsssnoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)


modelnetherlands_completenoother$timestamp.info$timestamps <- netherlands$year 
modelnetherlands_completenoother$timestamp.info$regular.timestamps <-netherlands$year 


###Country graphs

png("netherlands.png", width=17,height=17,units="cm", res=600)
plot(modelnetherlands_complete04,  ylim=c(0,1), scale="mean",
     #main="Netherlands",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

###Norway


norways<-read.csv2(file="Norwayshort.csv", fileEncoding = "ISO-8859-1")


norwaynew <- data.frame(year=norways$year[1]:norways$year[nrow(norways)]) 

norway <- full_join(norways,norwaynew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

norway$event <- with(norway, ifelse(is.na(event), 0, event))

norway$eventnoaftt<-with(norway, ifelse(year == 1871 | year == 1872 | year == 1893 & government == "Emil Stang" | year == 1894, NA, noncevent))

norway$noother<-with(norway, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


#SdPrior = 0.4

norwayss04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-5, 1))



modelnorway_complete04<- bsts(norway$noncevent,
                            state.specification=norwayss04,
                            niter=10000,
                            family="logit",
                            seed=123456,
)


modelnorway_complete04$timestamp.info$timestamps <- norway$year 
modelnorway_complete04$timestamp.info$regular.timestamps <-norway$year 

#With SdPrior = 0.3

norway03<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.3, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-5, 1)) 



modelnorway_complete03<- bsts(norway$noncevent,
                             state.specification=norway03,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelnorway_complete03$timestamp.info$timestamps <- norway$year 
modelnorway_complete03$timestamp.info$regular.timestamps <-norway$year 

#With SdPrior = 0.5

norway05<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.5, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-5, 1)) 



modelnorway_complete05<- bsts(norway$noncevent,
                             state.specification=norway05,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelnorway_complete05$timestamp.info$timestamps <- norway$year 
modelnorway_complete05$timestamp.info$regular.timestamps <-norway$year 

#Hosstart
norwayssJ<-AddLocalLevel(list(),  
                         sigma.prior=SdPrior(sigma.guess=.4, 
                                             sample.size=10,
                                             upper.limit=Inf),
                         initial.state.prior=NormalPrior(-5, 1))



modelnorway_completeJ<- bsts(norway$noncevent,
                             state.specification=norwayssJ,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelnorway_completeJ$timestamp.info$timestamps <- norway$year 
modelnorway_completeJ$timestamp.info$regular.timestamps <-norway$year 


#noaftt

norwayssnoaftt<-AddLocalLevel(list(),  
                              sigma.prior=SdPrior(sigma.guess=.4, 
                                                  sample.size=10,
                                                  upper.limit=Inf),
                              initial.state.prior=NormalPrior(-5, 1))



modelnorway_completenoaftt<- bsts(norway$eventnoaftt,
                                  state.specification=norwayssnoaftt,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelnorway_completenoaftt$timestamp.info$timestamps <- norway$year 
modelnorway_completenoaftt$timestamp.info$regular.timestamps <-norway$year 

#noother

norwayssnoother<-AddLocalLevel(list(),  
                              sigma.prior=SdPrior(sigma.guess=.4, 
                                                  sample.size=10,
                                                  upper.limit=Inf),
                              initial.state.prior=NormalPrior(-5, 1))



modelnorway_completenoother<- bsts(norway$noother,
                                  state.specification=norwayssnoother,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelnorway_completenoother$timestamp.info$timestamps <- norway$year 
modelnorway_completenoother$timestamp.info$regular.timestamps <-norway$year 


###Country graphs


png("norway.png", width=17,height=17,units="cm", res=600)
plot(modelnorway_complete04,  ylim=c(0,1), scale="mean",
     # main="Norway", 
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()


###Sweden


swedens<-read.csv2(file="Swedenshort.csv", fileEncoding = "ISO-8859-1")



swedennew <- data.frame(year=swedens$year[1]:swedens$year[nrow(swedens)]) 

sweden <- full_join(swedens,swedennew, by="year") %>% arrange(year) %>% mutate( event = ifelse(is.na(event), 0, event))

sweden$event <- with(sweden, ifelse(is.na(event), 0, event))

sweden$eventnoaftt <-with(sweden, ifelse(year == 1891 & nonpevent == 1, NA, noncevent))

sweden$noother<-with(sweden, ifelse(noncevent == 1 & pevent == 0, NA, pevent))


#SdPrior = 0.3

sweden03<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.3, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(2, 1)) 



modelsweden_complete03<- bsts(sweden$noncevent,
                             state.specification=sweden03,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelsweden_complete03$timestamp.info$timestamps <- sweden$year 
modelsweden_complete03$timestamp.info$regular.timestamps <-sweden$year 

#SdPrior = 0.4

sweden04<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(2, 1)) 



modelsweden_complete04<- bsts(sweden$noncevent,
                             state.specification=sweden04,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelsweden_complete04$timestamp.info$timestamps <- sweden$year 
modelsweden_complete04$timestamp.info$regular.timestamps <-sweden$year 

#SdPrior = 0.5

sweden05<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.5, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(2, 1)) 



modelsweden_complete05<- bsts(sweden$noncevent,
                             state.specification=sweden05,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelsweden_complete05$timestamp.info$timestamps <- sweden$year 
modelsweden_complete05$timestamp.info$regular.timestamps <-sweden$year 


#Hosstart
swedenJ<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-6, 1)) 



modelsweden_completeJ<- bsts(sweden$noncevent,
                             state.specification=swedenJ,
                             niter=10000,
                             family="logit",
                             seed=123456,
)


modelsweden_completeJ$timestamp.info$timestamps <- sweden$year 
modelsweden_completeJ$timestamp.info$regular.timestamps <-sweden$year 


#noaftt


swedenssnoaftt<-AddLocalLevel(list(),  
                              sigma.prior=SdPrior(sigma.guess=.4, 
                                                  sample.size=10,
                                                  upper.limit=Inf),
                              initial.state.prior=NormalPrior(2, 1)) 



modelsweden_completenoaftt<- bsts(sweden$eventnoaftt,
                                  state.specification=swedenssnoaftt,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelsweden_completenoaftt$timestamp.info$timestamps <- sweden$year 
modelsweden_completenoaftt$timestamp.info$regular.timestamps <-sweden$year 

###No other

swedenssnoother<-AddLocalLevel(list(),  
                              sigma.prior=SdPrior(sigma.guess=.4, 
                                                  sample.size=10,
                                                  upper.limit=Inf),
                              initial.state.prior=NormalPrior(2, 1)) 



modelsweden_completenoother<- bsts(sweden$noother,
                                  state.specification=swedenssnoother,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)


modelsweden_completenoother$timestamp.info$timestamps <- sweden$year 
modelsweden_completenoother$timestamp.info$regular.timestamps <-sweden$year 


###Country graphs

png("sweden.png", width=17,height=17,units="cm", res=600)
plot(modelsweden_complete04,  ylim=c(0,1), scale="mean",
     # main="Sweden",# sigma.guess=0.4",
     ylab="Parliamentarism", xlab="Year", show.actuals="TRUE")
dev.off()

####################
##Put it one graph##
####################

###Austria

meanpredsaustria03=modelaustria_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustria04=modelaustria_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustria05=modelaustria_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustriaconst=modelaustria_completeconst$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustriahosstart<-modelaustria_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustrianoaftt<-modelaustria_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsaustrianoother<-modelaustria_completenoother$state.contributions %>% as.data.frame() %>% colMeans()


austriamean03<-inv.logit(meanpredsaustria03)
austriamean04<-inv.logit(meanpredsaustria04)
austriamean05<-inv.logit(meanpredsaustria05)
austriameanconst<-inv.logit(meanpredsaustriaconst)
austriameanhosstart<-inv.logit(meanpredsaustriahosstart)
austriameannoaftt<-inv.logit(meanpredsaustrianoaftt)
austriameannoother<-inv.logit(meanpredsaustrianoother)


austria$country<-"Austria"

austria$mean03<-austriamean03
austria$mean03 <- ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$mean03)

austria$mean04<-austriamean04
austria$mean04 <- ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$mean04)

austria$mean05<-austriamean05
austria$mean05 <- ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$mean05)

austria$meanconst<-austriameanconst
austria$meanconst<-ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$meanconst)

austria$meanhosstart<-austriameanhosstart
austria$meanhosstart<-ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$meanhosstart)

austria$meannoaftt<-austriameannoaftt
austria$meannoaftt<-ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$meannoaftt)

austria$meannoother<-austriameannoother
austria$meannoother<-ifelse(austria$year > 1932 & austria$year < 1945, NA, austria$meannoother)



###Belgium

meanpredsbelgium03=modelbelgium_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgium04=modelbelgium_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgium05=modelbelgium_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgiumhosstart<-modelbelgium_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgiumnoaftt<-modelbelgium_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgiumnoaftt<-modelbelgium_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbelgiumnoother<-modelbelgium_completenoother$state.contributions %>% as.data.frame() %>% colMeans()



belgiummean03<-inv.logit(meanpredsbelgium03)
belgiummean04<-inv.logit(meanpredsbelgium04)
belgiummean05<-inv.logit(meanpredsbelgium05)
belgiummeanhosstart<-inv.logit(meanpredsbelgiumhosstart)
belgiummeannoaftt<-inv.logit(meanpredsbelgiumnoaftt)
belgiummeannoother<-inv.logit(meanpredsbelgiumnoother)


belgium$country<-"Belgium"

belgium$mean03<-belgiummean03

belgium$mean04<-belgiummean04

belgium$mean05<-belgiummean05

belgium$meanconst<-belgium$mean04

belgium$meanhosstart<-belgiummeanhosstart

belgium$meannoaftt<-belgiummeannoaftt

belgium$meannoother<-belgiummeannoother



###Britain

meanpredsbritain03=modelbritain_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbritain04=modelbritain_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbritain05=modelbritain_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbritainhosstart<-modelbritain_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbritainnoaftt<-modelbritain_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsbritainnoother<-modelbritain_completenoother$state.contributions %>% as.data.frame() %>% colMeans()



britainmean03<-inv.logit(meanpredsbritain03)
britainmean04<-inv.logit(meanpredsbritain04)
britainmean05<-inv.logit(meanpredsbritain05)
britainmeanhosstart<-inv.logit(meanpredsbritainhosstart)
britainmeannoaftt<-inv.logit(meanpredsbritainnoaftt)
britainmeannoother<-inv.logit(meanpredsbritainnoother)


britain$country<-"Britain"

britain$mean03<-britainmean03

britain$mean04<-britainmean04

britain$mean05<-britainmean05

britain$meanconst<-britain$mean04

britain$meanhosstart<-britainmeanhosstart

britain$meannoaftt<-britainmeannoaftt

britain$meannoother<-britainmeannoother


###Denmark

meanpredsdenmark03=modeldenmark_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsdenmark04=modeldenmark_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsdenmark05=modeldenmark_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsdenmarkhosstart=modeldenmark_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsdenmarknoaftt=modeldenmark_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsdenmarknoother=modeldenmark_completenoother$state.contributions %>% as.data.frame() %>% colMeans()



denmarkmean03<-inv.logit(meanpredsdenmark03)
denmarkmean04<-inv.logit(meanpredsdenmark04)
denmarkmean05<-inv.logit(meanpredsdenmark05)
denmarkmeanhosstart<-inv.logit(meanpredsdenmarkhosstart)
denmarkmeannoaftt<-inv.logit(meanpredsdenmarknoaftt)
denmarkmeannoother<-inv.logit(meanpredsdenmarknoother)


denmark$country<-"Denmark"

denmark$mean03<-denmarkmean03

denmark$mean04<-denmarkmean04

denmark$mean05<-denmarkmean05

denmark$meanconst<-denmarkmean04

denmark$meanhosstart<-denmarkmeanhosstart

denmark$meannoaftt<-denmarkmeannoaftt

denmark$meannoother<-denmarkmeannoother


###Finland

meanpredsfinland03=modelfinland_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfinland04=modelfinland_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfinland05=modelfinland_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfinlandhosstart=modelfinland_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfinlandnoother=modelfinland_completenoother$state.contributions %>% as.data.frame() %>% colMeans()



finlandmean03<-inv.logit(meanpredsfinland03)
finlandmean04<-inv.logit(meanpredsfinland04)
finlandmean05<-inv.logit(meanpredsfinland05)
finlandmeanhosstart<-inv.logit(meanpredsfinlandhosstart)
finlandmeannoother<-inv.logit(meanpredsfinlandnoother)


finland$country<-"Finland"

finland$mean03<-finlandmean03

finland$mean04<-finlandmean04

finland$mean05<-finlandmean05

finland$meanconst<-finlandmean05 #One major regime change in Finland

finland$meanhosstart<-finlandmeanhosstart

finland$meannoaftt<-finlandmean04 #No failed attempts in Finland

finland$meannoother<-finlandmeannoother


###France. Rerun models to recover state contributions
###SdPrior = 0.3

francess03<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.3, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete03<- bsts(france$noncevent,
                              state.specification=francess03,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modelfrance_complete03$timestamp.info$timestamps <- france$year 
modelfrance_complete03$timestamp.info$regular.timestamps <- france$year 

#With SdPrior=0.4

france04<-AddLocalLevel(list(),  
                        sigma.prior=SdPrior(sigma.guess=.4, 
                                            sample.size=10,
                                            upper.limit=Inf),
                        initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete04<- bsts(france$noncevent,
                              state.specification=france04,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modelfrance_complete04$timestamp.info$timestamps <- france$year 
modelfrance_complete04$timestamp.info$regular.timestamps <- france$year 

#With SdPrior=0.5

francess05<-AddLocalLevel(list(),  
                          sigma.prior=SdPrior(sigma.guess=.5, 
                                              sample.size=10,
                                              upper.limit=Inf),
                          initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_complete05<- bsts(france$noncevent,
                              state.specification=francess05,
                              niter=10000,
                              family="logit",
                              seed=123456,
)

modelfrance_complete05$timestamp.info$timestamps <- france$year 
modelfrance_complete05$timestamp.info$regular.timestamps <- france$year 

#France with SdPrior = 0.9

franceconst<-AddLocalLevel(list(),  
                           sigma.prior=SdPrior(sigma.guess=.9, 
                                               sample.size=10,
                                               upper.limit=Inf),
                           initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completeconst<- bsts(france$noncevent,
                                 state.specification=franceconst,
                                 niter=10000,
                                 family="logit",
                                 seed=123456,
)

modelfrance_completeconst$timestamp.info$timestamps <- france$year 
modelfrance_completeconst$timestamp.info$regular.timestamps <- france$year 

#Hosstart
franceJ<-AddLocalLevel(list(),  
                       sigma.prior=SdPrior(sigma.guess=.4, 
                                           sample.size=10,
                                           upper.limit=Inf),
                       initial.state.prior=NormalPrior(-7, 1)) 

modelfrance_completeJ<- bsts(france$noncevent,
                             state.specification=franceJ,
                             niter=10000,
                             family="logit",
                             seed=123456,
)

modelfrance_completeJ$timestamp.info$timestamps <- france$year 
modelfrance_completeJ$timestamp.info$regular.timestamps <- france$year 


#noaftt

francenoaftt<-AddLocalLevel(list(),  
                            sigma.prior=SdPrior(sigma.guess=.4, 
                                                sample.size=10,
                                                upper.limit=Inf),
                            initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completenoaftt<- bsts(france$eventnoaftt,
                                  state.specification=francenoaftt,
                                  niter=10000,
                                  family="logit",
                                  seed=123456,
)

modelfrance_completenoaftt$timestamp.info$timestamps <- france$year 
modelfrance_completenoaftt$timestamp.info$regular.timestamps <- france$year 


#noother

francenoother<-AddLocalLevel(list(),  
                             sigma.prior=SdPrior(sigma.guess=.4, 
                                                 sample.size=10,
                                                 upper.limit=Inf),
                             initial.state.prior=NormalPrior(-1, 1)) 

modelfrance_completenoother<- bsts(france$noother,
                                   state.specification=francenoother,
                                   niter=10000,
                                   family="logit",
                                   seed=123456,
)

modelfrance_completenoother$timestamp.info$timestamps <- france$year 
modelfrance_completenoother$timestamp.info$regular.timestamps <- france$year 


meanpredsfrance03=modelfrance_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfrance04=modelfrance_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfrance05=modelfrance_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfranceconst=modelfrance_completeconst$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfrancehosstart=modelfrance_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfrancenoaftt=modelfrance_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsfrancenoother=modelfrance_completenoother$state.contributions %>% as.data.frame() %>% colMeans()

francemean03<-inv.logit(meanpredsfrance03)
francemean04<-inv.logit(meanpredsfrance04)
francemean05<-inv.logit(meanpredsfrance05)
francemeanconst=inv.logit(meanpredsfranceconst)
francemeanhosstart=inv.logit(meanpredsfrancehosstart)
francemeannoaftt=inv.logit(meanpredsfrancenoaftt)
francemeannoother=inv.logit(meanpredsfrancenoother)

france$country<-"France"

france$mean03<-francemean03
france$mean03 <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$mean03)

france$mean04<-francemean04
france$mean04 <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$mean04)

france$mean05<-francemean05
france$mean05 <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$mean05)

france$meanconst<-francemeanconst
france$meanconst <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$meanconst)

france$meanhosstart<-francemeanhosstart
france$meanhosstart <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$meanhosstart)

france$meannoaftt<-francemeannoaftt
france$meannoaftt <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$meannoaftt)

france$meannoother<-francemeannoother
france$meannoother <- ifelse(france$year == 1940 & france$government == "Pétain" | france$year == 1940 & france$government == "Laval" | france$year > 1940 & france$year < 1945, NA, france$meannoother)

###Germany

meanpredsgermany03=modelgermany_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermany04=modelgermany_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermany05=modelgermany_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermanyconst=modelgermany_completeconst$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermanyhosstart=modelgermany_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermanynoaftt=modelgermany_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsgermanynoother=modelgermany_completenoother$state.contributions %>% as.data.frame() %>% colMeans()

germanymean03<-inv.logit(meanpredsgermany03)
germanymean04<-inv.logit(meanpredsgermany04)
germanymean05<-inv.logit(meanpredsgermany05)
germanymeanconst<-inv.logit(meanpredsgermanyconst)
germanymeanhosstart<-inv.logit(meanpredsgermanyhosstart)
germanymeannoaftt<-inv.logit(meanpredsgermanynoaftt)
germanymeannoother<-inv.logit(meanpredsgermanynoother)


germany$country<-"Germany"

germany$mean03<-germanymean03
germany$mean03 <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$mean03)

germany$mean04<-germanymean04
germany$mean04 <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$mean04)

germany$mean05<-germanymean05
germany$mean05 <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$mean05)

germany$meanconst<-germanymeanconst
germany$meanconst <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$meanconst)

germany$meanhosstart<-germanymeanhosstart
germany$meanhosstart <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$meanhosstart)

germany$meannoaftt<-germanymeannoaftt
germany$meannoaftt <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$meannoaftt)

germany$meannoother<-germanymeannoother
germany$meannoother <- ifelse(germany$year > 1933 & germany$year < 1949, NA, germany$meannoother)


###Italy

meanpredsitaly03=modelitaly_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsitaly04=modelitaly_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsitaly05=modelitaly_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsitalyhosstart=modelitaly_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsitalynoaftt=modelitaly_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsitalynoother=modelitaly_completenoother$state.contributions %>% as.data.frame() %>% colMeans()

italymean03<-inv.logit(meanpredsitaly03)
italymean04<-inv.logit(meanpredsitaly04)
italymean05<-inv.logit(meanpredsitaly05)
italymeanhosstart<-inv.logit(meanpredsitalyhosstart)
italymeannoaftt<-inv.logit(meanpredsitalynoaftt)
italymeannoother<-inv.logit(meanpredsitalynoother)


italy$country<-"Italy"

italy$mean03<-italymean03
italy$mean03 <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$mean03)

italy$mean04<-italymean04
italy$mean04 <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$mean04)

italy$mean05<-italymean05
italy$mean05 <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$mean05)

italy$meanconst<-italymean05
italy$meanconst <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$meanconst)

italy$meanhosstart<-italymeanhosstart
italy$meanhosstart <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$meanhosstart)

italy$meannoaftt<-italymeannoaftt
italy$meannoaftt <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$meannoaftt)

italy$meannoother<-italymeannoother
italy$meannoother <- ifelse(italy$year > 1922 & italy$year < 1945 | italy$year == 1945 & italy$government == "Bonomi", NA, italy$meannoother)




###Netherlands

meanpredsnetherlands03=modelnetherlands_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnetherlands04=modelnetherlands_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnetherlands05=modelnetherlands_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnetherlandshosstart=modelnetherlands_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnetherlandsnoother=modelnetherlands_completenoother$state.contributions %>% as.data.frame() %>% colMeans()



netherlandsmean03<-inv.logit(meanpredsnetherlands03)
netherlandsmean04<-inv.logit(meanpredsnetherlands04)
netherlandsmean05<-inv.logit(meanpredsnetherlands05)
netherlandsmeanhosstart<-inv.logit(meanpredsnetherlandshosstart)
netherlandsmeannoother<-inv.logit(meanpredsnetherlandsnoother)


netherlands$country<-"Netherlands"

netherlands$mean03<-netherlandsmean03

netherlands$mean04<-netherlandsmean04

netherlands$mean05<-netherlandsmean05

netherlands$meanconst<-netherlandsmean04 #no major regime changes in the Netherlands

netherlands$meanhosstart<-netherlandsmeanhosstart

netherlands$meannoaftt<-netherlandsmean04 #no failed attempts in the Netherlands

netherlands$meannoother<-netherlandsmeannoother #no failed attempts in the Netherlands



###Norway

meanpredsnorway03=modelnorway_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnorway04=modelnorway_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnorway05=modelnorway_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnorwayhosstart=modelnorway_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnorwaynoaftt=modelnorway_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsnorwaynoother=modelnorway_completenoother$state.contributions %>% as.data.frame() %>% colMeans()

norwaymean03<-inv.logit(meanpredsnorway03)
norwaymean04<-inv.logit(meanpredsnorway04)
norwaymean05<-inv.logit(meanpredsnorway05)
norwaymeanhosstart<-inv.logit(meanpredsnorwayhosstart)
norwaymeannoaftt<-inv.logit(meanpredsnorwaynoaftt)
norwaymeannoother<-inv.logit(meanpredsnorwaynoother)

norway$country<-"Norway"

norway$mean03<-norwaymean03

norway$mean04<-norwaymean04

norway$mean05<-norwaymean05

norway$meanconst<-norwaymean05

norway$meanhosstart<-norwaymeanhosstart

norway$meannoaftt<-norwaymeannoaftt

norway$meannoother<-norwaymeannoother



###Sweden

meanpredssweden03=modelsweden_complete03$state.contributions %>% as.data.frame() %>% colMeans()
meanpredssweden04=modelsweden_complete04$state.contributions %>% as.data.frame() %>% colMeans()
meanpredssweden05=modelsweden_complete05$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsswedenhosstart=modelsweden_completeJ$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsswedennoaftt=modelsweden_completenoaftt$state.contributions %>% as.data.frame() %>% colMeans()
meanpredsswedennoother=modelsweden_completenoother$state.contributions %>% as.data.frame() %>% colMeans()


swedenmean03<-inv.logit(meanpredssweden03)
swedenmean04<-inv.logit(meanpredssweden04)
swedenmean05<-inv.logit(meanpredssweden05)
swedenmeanhosstart<-inv.logit(meanpredsswedenhosstart)
swedenmeannoaftt<-inv.logit(meanpredsswedennoaftt)
swedenmeannoother<-inv.logit(meanpredsswedennoother)


sweden$country<-"Sweden"

sweden$mean03<-swedenmean03

sweden$mean04<-swedenmean04

sweden$mean05<-swedenmean05

sweden$meanconst<-swedenmean05

sweden$meanhosstart<-swedenmeanhosstart

sweden$meannoaftt<-swedenmeannoaftt

sweden$meannoother<-swedenmeannoother



Parldev<-rbind.fill(list(austria, belgium, britain, denmark, finland, france, germany, italy, netherlands, norway, sweden))
#Parldev<- Parldev %>%
 # mutate(country2=country)

syncdata2<- Parldev %>%
  mutate(country2=country)

###Main graph

syncdata2 %>%
  ggplot( aes(y=mean04, x=year, group=country, color=country)) +
  geom_line( data=syncdata2 %>% dplyr::select(-country), aes(group=country2), color="grey", size=0.5, alpha=0.5) +
  geom_line( aes(color=country), color="#150f20", size=1.0 )+
  scale_color_viridis(discrete = TRUE) +
  #theme_ipsum() +
  theme(
    legend.position="none",
    plot.title = element_text(size=14),
    panel.grid = element_blank()) +
  ggtitle("") + 
  ylab("Parliamentarism")+
  xlab("Year")+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank())+ scale_colour_discrete(name="Country")+
  theme(plot.title = element_text(hjust = 0.5))+
  facet_wrap(~country) 

ggsave("trajectoriescomparednotrend0.4.pdf", width = 24, height = 16, units = "cm")
ggsave("trajectoriescomparednotrend0.4.png", width = 24, height = 16, units = "cm", dpi = 600)

remove(syncdata2)


#Adjusting for major regime changes ("const")

trajectoriescomparednotrendbase0.4<-ggplot(data=Parldev, aes(year)) + geom_line(aes(y=mean04, colour="mean04")) +
  geom_line(aes(y=meanconst, colour="meanconst")) + 
  facet_wrap(~country)+
  labs(title = "", x = "Year", y = "", color = "") + ylim(0,1)+
  scale_color_manual(labels = c("Parliamentarism", "Parliamentarism\nadjusted for major\nregime changes"), values = c("black", "lightgray")) +
  theme(legend.position =c(0.88, 0.15))+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))+xlab("Year")+ylab("")

trajectoriescomparednotrendbase0.4=trajectoriescomparednotrendbase0.4 + theme(legend.position =c (0.88, 0.15))


ggsave("trajectoriescomparednotrendbase0.4.pdf", width = 24, height = 16, units = "cm")
ggsave("trajectoriescomparednotrendbase0.4.png", width = 24, height = 16, units = "cm", dpi=600)


#Comparing main results with trajectories excluding "other" resignations (noother)

prulenoother<-ggplot(data=Parldev, aes(year)) + geom_line(aes(y=mean04, colour="mean04")) +
  geom_line(aes(y=meannoother, colour="meannoother")) + 
  facet_wrap(~country)+
  labs(title = "", x = "Year", y = "", color = "") + ylim(0,1)+
  scale_color_manual(labels = c("Parliamentarism", "Parliamentarism\nexcluding 'other'\nresignations"), values = c("black", "lightgray")) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))+xlab("Year")+ylab("")

prulenoother=prulenoother + theme(legend.position =c (0.88, 0.15))

ggsave("prule+noother.pdf", width = 24, height = 16, units = "cm")
ggsave("prule+noother.png", width = 24, height = 16, units = "cm", dpi=600)


#Comparing main results with trajectories that exclude failed attempts and that takes heads of state to dominate at the establishment of parliaments

prulehosstartnoaftt<-ggplot(data=Parldev, aes(year)) + geom_line(aes(y=mean04, colour="mean04")) +
    geom_line(aes(y=meanhosstart, colour="meanhosstart")) + 
  geom_line(aes(y=meannoaftt, colour="meannoaftt")) + 
  facet_wrap(~country)+
  labs(title = "", x = "Year", y = "", color = "") + ylim(0,1)+
  scale_color_manual(labels = c("Parliamentarism", "Parliamentarism\nwith initial head-of-\nstate dominance", "Parliamentarism\nwithout failed attempts"), values = c("black", "darkgray", "lightgray")) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))+xlab("Year")+ylab("")

prulehosstartnoaftt=prulehosstartnoaftt+theme(legend.position =c (0.88, 0.13))

ggsave("prule+hosstart+noaftt.pdf", width = 24, height = 16, units = "cm")
ggsave("prule+hosstart+noaftt.png", width = 24, height = 16, units = "cm", dpi=600)


###Comparing main graph with sigma = 0.3 and sigma = 0.5

sigmaadjustments<-ggplot(data=Parldev, aes(year)) + 
  geom_line(aes(y=mean03, colour="mean03")) +
  geom_line(aes(y=mean04, colour="mean04")) + 
  geom_line(aes(y=mean05, colour="mean05")) +
  facet_wrap(~country)+
  labs(title = "", x = "Year", y = "", color = "") + ylim(0,1)+
  scale_color_manual(labels = c("Parliamentarism, s=0.3", "Parliamentarism, s=0.4", "Parliamentarism, s=0.5"), 
                     values = c("darkgray", "black", "lightgray")) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))+xlab("Year")+ylab("")

sigmaadjustments=sigmaadjustments+theme(legend.position =c (0.885, 0.13))
  
ggsave("sigmaadjustments.pdf", width = 24, height = 16, units = "cm")
ggsave("sigmaadjustments.png", width = 24, height = 16, units = "cm", dpi=600)


####List of variables produced

### mean04=parliamentarism, main estimation
### mean03=parliamentarism, decreased uncertainty paramenter
### mean05=parliamentarism, increased uncertainty parameter
### meanconst = estimation adjusted for major constitutional reforms or regime changes
### meanhosstart = estimation assuming all countries were personal monarchies at the start,
# i.e., that only the head of state is expected to terminate the government
### meannoaftt=estimation that disregards failed termination attempts (no other actor fails to terminate)
### meannoother=estimation that disregards "other" events, i.e., those in which neither parliament
# nor the head of state is involved.

#Write out data
library(openxlsx)
write.xlsx(Parldev, file="Parldev.xlsx")

#Comparing to V-dem confidence index

setwd("C:/Users/simda81/OneDrive - Linköpings universitet/Parliamentary Development/Datasets/2020")


vdem<-readRDS("V-Dem-CY-Full+Others-v10.rds")

vdem <- dplyr::select(vdem, country_name, year, v2x_ex_confidence)

vdem<-vdem %>% dplyr::mutate(country_name = case_when(
  country_name == "United Kingdom"~"Britain",
  country_name == "Piedmont-Sardinia"~"Italy",
  TRUE~country_name
))

vdem <- vdem %>%  dplyr::mutate(country = country_name)

Parldev<-Parldev[-c(1567),]


Parldev<-left_join(Parldev, vdem, join_by(country, year))

Parldev<-Parldev[-c(1568),]


remove(vdem)

setwd("C:/Users/simda81/OneDrive - Linköpings universitet/Parliamentary Development/Simon avhandling/BJPS/Replication files and data")

pruleagainstvdemOLS<-ggplot(data=Parldev, aes(year)) + geom_line(aes(y=mean04, colour="mean04")) + 
  geom_line(aes(y=v2x_ex_confidence, colour="v2x_ex_confidence")) + facet_wrap(~country)+
  labs(title = "", x = "Year", y = "", color = "") + ylim(0,1)+
  scale_color_manual(labels = c("Parliamentarism", "V-Dem's Confidence\nIndex\n(v2x_ex_confidence)"), values = c("black", "lightgray")) +
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "black"))+xlab("Year")+ylab("")

pruleagainstvdemOLS=pruleagainstvdemOLS+theme(legend.position =c(0.88, 0.13))

ggsave("pruleagainstvdemOLS.pdf", width = 24, height = 16, units = "cm")
ggsave("pruleagainstvdemOLS.png", width = 24, height = 16, units = "cm", dpi=600)


###Correlations

library(Hmisc)

rcorr(Parldev$mean04, Parldev$v2x_ex_confidence, type="pearson") #0.66, p=0, n=2296
rcorr(Parldev$mean04, Parldev$meannoother, type="pearson") #0.94, p=0, n=2415
rcorr(Parldev$mean04, Parldev$meannoaftt, type="pearson") #0.98, p=0, n=2415
rcorr(Parldev$mean04, Parldev$meanconst, type="pearson") #0.99, p=0, n=2415
rcorr(Parldev$mean04, Parldev$meanhosstart, type="pearson") #0.95, p=0, n=2415


